#
# step10 Cluster the cells
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)



#
()
#./seurat-4.1.0/R/generics.R:179:FindNeighbors <- function(object, ...) {

对于给定的数据集，计算 k 个最近的邻居。
通过计算每个细胞的临近的重叠情况(Jaccard 指数)，以及它的 k 个最近的邻居，
可以 可选的通过 compute.SNN 构建一个共享最近邻图，


#' (Shared) Nearest-neighbor graph construction # (共享)最近邻图构建
#'
#' Computes the \code{k.param} nearest neighbors for a given dataset. Can also
#' optionally (via \code{compute.SNN}), construct a shared nearest neighbor
#' graph by calculating the neighborhood overlap (Jaccard index) between every
#' cell and its \code{k.param} nearest neighbors.
#'
#' @param object An object
#' @param ... Arguments passed to other methods
#'
#' @return This function can either return a \code{\link{Neighbor}} object
#' with the KNN information or a list of \code{\link{Graph}} objects with
#' the KNN and SNN depending on the settings of \code{return.neighbor} and
#' \code{compute.SNN}. When running on a \code{\link{Seurat}} object, this
#' returns the \code{\link{Seurat}} object with the Graphs or Neighbor objects
#' stored in their respective slots. Names of the Graph or Neighbor object can
#' be found with \code{\link{Graphs}} or \code{\link{Neighbors}}.
#'
#' @examples
#' data("pbmc_small")
#' pbmc_small
#' # Compute an SNN on the gene expression level #使用高变基因
#' pbmc_small <- FindNeighbors(pbmc_small, features = VariableFeatures(object = pbmc_small))
#'
#' # More commonly, we build the SNN on a dimensionally reduced form of the data
#' # such as the first 10 principle components. #使用降维后的前10个PC
#'
#' pbmc_small <- FindNeighbors(pbmc_small, reduction = "pca", dims = 1:10)
#'
#' @rdname FindNeighbors
#' @export FindNeighbors
#'
FindNeighbors <- function(object, ...) {
  UseMethod(generic = 'FindNeighbors', object = object)
}






()
#./seurat-4.1.0/R/clustering.R:729:FindNeighbors.Seurat <- function(

#' @param assay Assay to use in construction of (S)NN; used only when \code{dims}
#' is \code{NULL} #用于构建SNN的 assay，仅在 dims为 NULL 时有效;
#' @param features Features to use as input for building the (S)NN; used only when
#' \code{dims} is \code{NULL} #用于构建 SNN的基因，仅在 dims为NULL时有效。
#' @param reduction Reduction to use as input for building the (S)NN #构建SNN时，输入哪个降维，默认是pca
#' @param dims Dimensions of reduction to use as input # 输入降维的哪些维度，默认 1:10
#' @param do.plot Plot SNN graph on tSNE coordinates #把 SNN 图输出到 tSNE 坐标上
#' @param graph.name Optional naming parameter for stored (S)NN graph #保存SNN图的 可选名字参数
#' (or Neighbor object, if return.neighbor = TRUE). Default is assay.name_(s)nn. #(如果 return.neighbor = TRUE，则是 Neighbor 对象 )
#' #默认是 assay.name_(s)nn
#' To store both the neighbor graph and the shared nearest neighbor (SNN) graph, #为了同时保存 neighbor 图和 SNN 图
#' you must supply a vector containing two names to the \code{graph.name} #则必须提供一个向量，包含graph.name的2个名字
#' parameter. 
#' The first element in the vector will be used to store the nearest # 第一个元素用于保存 NN 图，第二个保存SNN图。
#' neighbor (NN) graph, and the second element used to store the SNN graph. 
#' If only one name is supplied, only the NN graph is stored. #如果只提供一个元素，则只保存 NN 图。
#'
#' @importFrom igraph graph.adjacency plot.igraph E # igraph 包擅长做什么？
#'
#' @rdname FindNeighbors
#' @export
#' @concept clustering
#' @method FindNeighbors Seurat
#'
FindNeighbors.Seurat <- function(
  object,
  reduction = "pca",
  dims = 1:10,
  assay = NULL,
  features = NULL,
  k.param = 20,  #影响什么？怎么设置？

  return.neighbor = FALSE,
  compute.SNN = !return.neighbor, #这两个是互斥的

  prune.SNN = 1/15, #剪枝？
  nn.method = "annoy", #最近邻算法
  n.trees = 50, #
  annoy.metric = "euclidean", #欧拉距离

  nn.eps = 0, 
  verbose = TRUE,
  force.recalc = FALSE,

  do.plot = FALSE, #默认不画图
  graph.name = NULL,
  l2.norm = FALSE, #L2 标准化
  cache.index = FALSE, #是否缓存 index
  ...
) {
  #(A1) 检查其他参数 //todo
  CheckDots(...)

  #(A2) 如果 dims 非空, 默认走这里。
  # pbmc <- FindNeighbors(pbmc, dims = 1:10) # 指定 dims 时，就是走的这里。
  if (!is.null(x = dims)) {
  	#(B1) 取出实验 DefaultAssay(pbmc_small[["pca"]]) #[1] "RNA"
    assay <- DefaultAssay(object = object[[reduction]])
    #(B2) 取出细胞的PC坐标
    # > class(Embeddings(pbmc_small[["pca"]]))
	# [1] "matrix" "array"
    data.use <- Embeddings(object = object[[reduction]])

    #(B3) 如果 输入参数dims最大值 大于 PC矩阵的列数，则报错
    if (max(dims) > ncol(x = data.use)) {
      stop("More dimensions specified in dims than have been computed")
    }

    #(B4) 取子集：PC的某些列，行为细胞
    data.use <- data.use[, dims]

    #(B5) 丢给了 FindNeighbors.default?
    neighbor.graphs <- FindNeighbors(
      object = data.use, #这是一个 "matrix"，行为cell，列为 PC
      k.param = k.param,
      compute.SNN = compute.SNN,
      prune.SNN = prune.SNN,
      nn.method = nn.method,
      n.trees = n.trees,
      annoy.metric = annoy.metric, #默认euclidean
      nn.eps = nn.eps,
      verbose = verbose,
      force.recalc = force.recalc,
      l2.norm = l2.norm, #默认F
      return.neighbor = return.neighbor,
      cache.index = cache.index,
      ...
    )

  # 如果 dims 是空的，默认不走这里
  } else {
  	# 获取默认实验名
    assay <- assay %||% DefaultAssay(object = object)
    # 获取 实验对象
    data.use <- GetAssay(object = object, assay = assay)
    
    # 丢给了 FindNeighbors.Assay()
    neighbor.graphs <- FindNeighbors(
      object = data.use,
      features = features,
      k.param = k.param,
      compute.SNN = compute.SNN,
      prune.SNN = prune.SNN,
      nn.method = nn.method,
      n.trees = n.trees,
      annoy.metric = annoy.metric,
      nn.eps = nn.eps,
      verbose = verbose,
      force.recalc = force.recalc,
      l2.norm = l2.norm,
      return.neighbor = return.neighbor,
      cache.index = cache.index,
      ...
    )
  }

  #(A3) 如果 neighbor.graphs 长度为1，则把它放到 list(nn=) 中
  # 默认长度是2，默认不走这里
  if (length(x = neighbor.graphs) == 1) {
    neighbor.graphs <- list(nn = neighbor.graphs)
  }


  #(A4) graph.name 如果没有提供，则设置默认值
  graph.name <- graph.name %||%
    # 如果有 return.neighbor=T，则使用.号连接
    if (return.neighbor) {
      paste0(assay, ".", names(x = neighbor.graphs))
    # 否则(默认)，则使用_连接
    } else {
      paste0(assay, "_", names(x = neighbor.graphs))
    }

  #(A5) 如果 graph.name 长度为1，则提示: 仅提供一个名字，则只保存NN(而不保存SNN)
  if (length(x = graph.name) == 1) {
    message("Only one graph name supplied, storing nearest-neighbor graph only")
  }

  #(A6) 遍历 graph.name
  for (ii in 1:length(x = graph.name)) {
  	# 如果 其中的元素 是 Graph 类，则设置其默认assay为 输入参数 assay
    if (inherits(x = neighbor.graphs[[ii]], what = "Graph")) {
      DefaultAssay(object = neighbor.graphs[[ii]]) <- assay
    }
    # 把该对象保存到 Seurat 对象的 graph.name 中 //写入 Seurat 对象
    # > names(pbmc@graphs)
	# [1] "RNA_nn"  "RNA_snn"
	#> class(pbmc@graphs$RNA_nn)
	#[1] "Graph"
    object[[graph.name[[ii]]]] <- neighbor.graphs[[ii]]
  }


  # (A7) 默认是不画图
  if (do.plot) {
  	# 如果没有 tsne，则提醒先运行 RunTSNE()
    if (!"tsne" %in% names(x = object@reductions)) {
      warning("Please compute a tSNE for SNN visualization. See RunTSNE().")
    # 否则，可以画图
    } else {
      # (B1) tSNE 的行数(cid) 如果 不等于 Seurat 的列数，则警告
      if (nrow(x = Embeddings(object = object[["tsne"]])) != ncol(x = object)) {
        warning("Please compute a tSNE for SNN visualization. See RunTSNE().")
      } else {
      	# 画网络图，根据 SNN 每个细胞的连接情况，默认使用 tSNE 坐标。
      	#(C1) 这个图见本文 3.2
        net <- graph.adjacency(
          adjmatrix = as.matrix(x = neighbor.graphs[[2]]),
          mode = "undirected",
          weighted = TRUE,
          diag = FALSE
        )
        #(C2) 这个图见本文 3.2
        plot.igraph(
          x = net,
          layout = as.matrix(x = Embeddings(object = object[["tsne"]])),
          edge.width = E(graph = net)$weight,
          vertex.label = NA,
          vertex.size = 0
        )
      }
    }
  }

  # (A8) 记录日志
  object <- LogSeuratCommand(object = object)
  return(object)
}








()FindNeighbors.default()
#./seurat-4.1.0/R/clustering.R:516:FindNeighbors.default <- function(



#' @param query Matrix of data to query against object. If missing, defaults to
#' object.
#' @param distance.matrix Boolean value of whether the provided matrix is a
#' distance matrix; note, for objects of class \code{dist}, this parameter will
#' be set automatically #是不是距离矩阵？如果是 dist 类，则该参数会自动设置
#' @param k.param Defines k for the k-nearest neighbor algorithm #定义 k 最近邻算法的 k
#' @param return.neighbor Return result as \code{\link{Neighbor}} object. Not
#' used with distance matrix input. #返回结果作为 Neighbor 对象。 对于输入 距离矩阵的 该参数不适用。
#' @param compute.SNN also compute the shared nearest neighbor graph #计算SNN
#' @param prune.SNN Sets the cutoff for acceptable Jaccard index when #构建SNN时接受的 Jaccard 索引的截断点。
#' computing the neighborhood overlap for the SNN construction. 
#' Any edges with 
#' values less than or equal to this will be set to 0 and removed from the SNN
#' graph. #<= 该值的边将会从 SNN 中去掉
#' Essentially sets the stringency of pruning (0 --- no pruning, 1 ---
#' prune everything). #设置剪枝的严谨程度， 0-不修剪，1-修剪
#' @param nn.method Method for nearest neighbor finding. Options include: rann,
#' annoy #最近邻算法的种类，可选: rann, annoy
#' @param annoy.metric Distance metric for annoy. Options include: euclidean,
#' cosine, manhattan, and hamming #几种距离的计算，可选 euclidean, cosine, manhattan, and hamming
#' @param n.trees More trees gives higher precision when using annoy approximate
#' nearest neighbor search #做最近邻搜索时，更多的树，就有更高的分辨率。
#' @param nn.eps Error bound when performing nearest neighbor seach using RANN;
#' default of 0.0 implies exact nearest neighbor search #使用RANN做最近邻搜索时的错误极限。
#' @param verbose Whether or not to print output to the console
#' @param force.recalc Force recalculation of (S)NN. #强制重新计算 SNN
#' @param l2.norm Take L2Norm of the data #计算数据的 L2 标准化
#' @param cache.index Include cached index in returned Neighbor object
#' (only relevant if return.neighbor = TRUE) #
#' @param index Precomputed index. Useful if querying new data against existing
#' index to avoid recomputing. 提起计算索引。当查询新数据 against 一直索引时避免重新计算，特别有效。
#'
#' @importFrom RANN nn2
#' @importFrom methods as
#'
#' @rdname FindNeighbors
#' @export
#' @concept clustering
#' @method FindNeighbors default
#'
FindNeighbors.default <- function(
  object,
  query = NULL,
  distance.matrix = FALSE,
  k.param = 20,

  return.neighbor = FALSE,
  compute.SNN = !return.neighbor,

  prune.SNN = 1/15,
  nn.method = "annoy",
  n.trees = 50,
  annoy.metric = "euclidean",
  nn.eps = 0,
  verbose = TRUE,
  force.recalc = FALSE,
  l2.norm = FALSE,
  cache.index = FALSE,
  index = NULL,
  ...
) {
  #(A1) 检查参数
  CheckDots(...)
  #(A2) 如果 obj 的dim()为空，则警告：强转为矩阵
  if (is.null(x = dim(x = object))) {
    warning(
      "Object should have two dimensions, attempting to coerce to matrix",
      call. = FALSE
    )
    # 强转为矩阵
    object <- as.matrix(x = object)
  }

  #(A3) 如果行名为空，则报错
  # 默认 object 是一个 "matrix"，行为cell，列为 PC
  if (is.null(rownames(x = object))) {
    stop("Please provide rownames (cell names) with the input object")
  }

  #(A4) 细胞个数 
  n.cells <- nrow(x = object)
  #如果细胞数 小于 k，则警告，并设置k为细胞数-1
  if (n.cells < k.param) {
    warning(
      "k.param set larger than number of cells. Setting k.param to number of cells - 1.",
      call. = FALSE
    )
    k.param <- n.cells - 1
  }

  #(A5) 如果 l2.norm =T，默认是F
  if (l2.norm) {
  	# 对矩阵进行 L2 标准化
    object <- L2Norm(mat = object)
    # `%iff%` 如果x非空则返回y，如果空则使用x。
    # 默认 query 是空，则保持空
    query <- query %iff% L2Norm(mat = query)
  }

  #(A6) 如果空，则使用矩阵 obj
  query <- query %||% object


  #(A7) 找到每个细胞的 k最近邻
  # find the k-nearest neighbors for each single cell
  #是否距离矩阵，默认否
  if (!distance.matrix) { 
    if (verbose) { #默认T
      if (return.neighbor) { #默认F，走else
        message("Computing nearest neighbors")
      } else {
        message("Computing nearest neighbor graph")
      }
    }

    #(B1)
    ##################
    # 干活的!
    ##################
    nn.ranked <- NNHelper(
      data = object,
      query = query,
      k = k.param,
      method = nn.method,
      n.trees = n.trees,
      searchtype = "standard",
      eps = nn.eps,
      metric = annoy.metric,
      cache.index = cache.index, #默认不保存索引
      index = index
    )
 	
 	# (B2) 默认F, 不走这里
    if (return.neighbor) {
      if (compute.SNN) {
        warning("The SNN graph is not computed if return.neighbor is TRUE.", call. = FALSE)
      }
      return(nn.ranked)
    }

    # (B3) 获取最近邻细胞的行号构成的矩阵
    nn.ranked <- Indices(object = nn.ranked)

  #是否距离矩阵，默认否。如果用户设置T才走这里
  } else {
  	# 提示：基于提供的距离矩阵，构建SNN
    if (verbose) {
      message("Building SNN based on a provided distance matrix")
    }
    # 空矩阵，n行(细胞数) x k列(k个最近邻的行号)
    knn.mat <- matrix(data = 0, ncol = k.param, nrow = n.cells)
    # 另一个空矩阵 ?
    knd.mat <- knn.mat
    # 遍历细胞
    for (i in 1:n.cells) {
      # 排序，最小的下标排在最前面，取前k个
      knn.mat[i, ] <- order(object[i, ])[1:k.param]
      # 把最小的k个距离下标，换成距离的值
      knd.mat[i, ] <- object[i, knn.mat[i, ]] #这个有啥用？后面没再用了 //todo
    }
    # 距离的值的前k列，作为每个点最近邻的k个点的下标
    nn.ranked <- knn.mat[, 1:k.param]
  }
  

  #(A8) 把最近邻点的下标转为图
  # convert nn.ranked into a Graph
  # 转置后变为 一列是一个点的最近邻的k个下标，强转为数组，df按列合并为一个数组。
  # 每k个就是一个点的最近邻点。
  j <- as.numeric(x = t(x = nn.ranked))
  # > ((1:10)-1) %/% 3 +1 #把1:10分组，每3个一组
  #[1] 1 1 1 2 2 2 3 3 3 4
  i <- ((1:length(x = j)) - 1) %/% k.param + 1


  #(A9) 定义和矩阵同行列的稀疏矩阵: 第i行(细胞等同的行)，第j列(最近邻的k个点)
  # 强转为 Graph 结构: 是 dgCMatrix 的子类，包含一个 slot: g@assay.used
  nn.matrix <- as(object = sparseMatrix(i = i, j = j, x = 1, dims = c(nrow(x = object), nrow(x = object))), 
  	Class = "Graph")
  # 增加行名、列名
  rownames(x = nn.matrix) <- rownames(x = object)
  colnames(x = nn.matrix) <- rownames(x = object)
  # 近邻图就是这个稀疏矩阵，每行是cell，其最近邻的点用1表示
  neighbor.graphs <- list(nn = nn.matrix)


  #(A10) 是否计算 SNN，默认是T
  if (compute.SNN) {
    if (verbose) {
      message("Computing SNN")
    }

    ####################
    # 计算SNN的主力函数，返回矩阵是雅各布
    ####################
    snn.matrix <- ComputeSNN(
      nn_ranked = nn.ranked, #每个细胞一行，一行是最近邻的k个点的下标
      prune = prune.SNN #剪枝，默认是 1/15
    )
    # 添加行名列名
    rownames(x = snn.matrix) <- rownames(x = object)
    colnames(x = snn.matrix) <- rownames(x = object)
    # 强转为 Graph 类
    snn.matrix <- as.Graph(x = snn.matrix)
    # 为list添加一个新的键值对
    neighbor.graphs[["snn"]] <- snn.matrix
  }

  return(neighbor.graphs)
}







() NNHelper()
#./seurat-4.1.0/R/clustering.R:1599:NNHelper <- function(data, query = data, k, method, cache.index = FALSE, ...) {

内部帮助类，实现 多种近邻查找方法的分发，目前就2种方法

# Internal helper function to dispatch to various neighbor finding methods
#
# @param data Input data
# @param query Data to query against data
# @param k Number of nearest neighbors to compute
# @param method Nearest neighbor method to use: "rann", "annoy"
# @param cache.index Store algorithm index with results for reuse
# @param ... additional parameters to specific neighbor finding method
#
#' @importFrom methods new
#' @importClassesFrom SeuratObject Neighbor
#
NNHelper <- function(data, query = data, k, method, cache.index = FALSE, ...) {
  # 获取本层的环境变量？ 见本文 3.7, 3.8
  args <- as.list(x = sys.frame(which = sys.nframe()))
  # 获取 ... 参数
  args <- c(args, list(...))

  # 根据 method 字符串，决定使用哪个函数？
  results <- (
    switch(
      EXPR = method,
      # 不是默认方法
      "rann" = {
        args <- args[intersect(x = names(x = args), y = names(x = formals(fun = nn2)))]
        do.call(what = 'nn2', args = args)
      },
      # 默认是"annoy"
      "annoy" = {
      	# 获取参数列表的名字(names())，和 求 NN的函数的参数名字求交集，然后获取这些参数及其值，见本文 3.9
        args <- args[intersect(x = names(x = args), y = names(x = formals(fun = AnnoyNN)))]
        do.call(what = 'AnnoyNN', args = args) #调用 what 函数，使用args 参数列表
      },
      # 求最近邻的方法，目前只能二选一，否则报错
      stop("Invalid method. Please choose one of 'rann', 'annoy'")
    )
  )

  # 新建 Neighbor类的对象，填充 上一步的结果
  n.ob <- new(
    Class = 'Neighbor',
    nn.idx = results$nn.idx,
    nn.dist = results$nn.dists,
    alg.info = results$alg.info %||% list(),
    cell.names = rownames(x = query)
  )

  # 如果 cache.index==T 且 results$idx 非空，则保存该数据：最近邻建立的树文件，作为索引
  if (isTRUE(x = cache.index) && !is.null(x = results$idx)) {
    slot(object = n.ob, name = "alg.idx") <- results$idx
  }
  
  return(n.ob)
}











() L2Norm()
# ./seurat-4.1.0/R/utilities.R:1692:L2Norm <- function(mat, MARGIN = 1){

详情见本文 3.5.







() AnnoyNN()
# ./seurat-4.1.0/R/clustering.R:858:AnnoyNN <- function(data,


# Run annoy
#
# @param data Data to build the index with
# @param query A set of data to be queried against data
# @param metric Distance metric; can be one of "euclidean", "cosine", "manhattan",
# "hamming"
# @param n.trees More trees gives higher precision when querying
# @param k Number of neighbors
# @param search.k During the query it will inspect up to search_k nodes which
# gives you a run-time tradeoff between better accuracy and speed.
# @param include.distance Include the corresponding distances
# @param index optional index object, will be recomputed if not provided
#
AnnoyNN <- function(data,
                    query = data,
                    metric = "euclidean",
                    n.trees = 50,
                    k,
                    search.k = -1,
                    include.distance = TRUE,
                    index = NULL
) {
  #(A1) 第一步构建索引 AnnoyBuildIndex()
  idx <- index %||% AnnoyBuildIndex(
    data = data,
    metric = metric,
    n.trees = n.trees)

  #(A2) 第二步 查找索引
  nn <- AnnoySearch(
    index = idx,
    query = query,
    k = k,
    search.k = search.k,
    include.distance = include.distance)

  #(A3) 加相关数据
  # 加上索引
  nn$idx <- idx
  # 求距离的方式，使用了几个PC
  nn$alg.info <- list(metric = metric, ndim = ncol(x = data))

  #(A4) 返回最近邻
  return(nn)
}











()AnnoyBuildIndex()

# Build the annoy index #这个函数内，作者的变量命名十分不讲究。
#
# @param data Data to build the index with
# @param metric Distance metric; can be one of "euclidean", "cosine", "manhattan",
# "hamming"
# @param n.trees More trees gives higher precision when querying
#
#' @importFrom RcppAnnoy AnnoyEuclidean AnnoyAngular AnnoyManhattan AnnoyHamming
#
AnnoyBuildIndex <- function(data, metric = "euclidean", n.trees = 50) {
  # (A1)这个 data 就是 cell行 x PC列
  f <- ncol(x = data)
  # (A2)选择要调用的C++代码的函数名，默认欧氏距离
  a <- switch(
    EXPR = metric,
    "euclidean" =  new(Class = RcppAnnoy::AnnoyEuclidean, f), #默认

    "cosine" = new(Class = RcppAnnoy::AnnoyAngular, f),
    "manhattan" = new(Class = RcppAnnoy::AnnoyManhattan, f),
    "hamming" = new(Class = RcppAnnoy::AnnoyHamming, f),
    stop ("Invalid metric")
  )

  # (A3)ii变量，逐行遍历细胞行号 1:nrow(data)
  for (ii in seq(nrow(x = data))) {
  	# 这个是啥方法？
    a$addItem(ii - 1, data[ii, ])
  }

  # (A4)构建树的个数，默认 50个
  a$build(n.trees)

  # 测试
  #> class(a)
  # [1] "Rcpp_AnnoyEuclidean"
  #> a
  # C++ object <0x55881a1c9bd0> of class 'AnnoyEuclidean' <0x5587fb478ef0>
  #> str(a)
  # Reference class 'Rcpp_AnnoyEuclidean' [package "RcppAnnoy"] with 0 fields
  # 看样子是一个 RC/R5 类，怎么看该类的内部保存的数据呢? 看搜索函数能不能提供线索
  return(a)
}










() RcppAnnoy::AnnoyEuclidean() //todo
'Annoy' 是一个C++库，计算 近似最近邻的。这里基于Rcpp包提供的R接口和Python包装类一致: https://github.com/spotify/annoy
https://cran.r-project.org/web/packages/RcppAnnoy/index.html

这里难度有点大，参考一下其他人的心得体会
https://www.jianshu.com/p/72313d70d9ab

该算法需要2步: 建立索引，搜索最近的点。
RcppAnnoy 包实现了这个过程。 不想看这个包的源码了，先放过。 //todo


> RcppAnnoy::AnnoyEuclidean
C++ class 'AnnoyEuclidean' <0x55d6b60e4000>
Constructors:
    AnnoyEuclidean(int)
        docstring : constructor with integer count

Fields: No public fields exposed by this class

Methods: 
     void addItem(int, Rcpp::NumericVector)  
           docstring : add item
     void build(int)  
           docstring : build an index
##############










() AnnoySearch()
./seurat-4.1.0/R/clustering.R:924:AnnoySearch <- function(index, query, k, search.k = -1, include.distance = TRUE) {


# Search an Annoy approximate nearest neighbor index #KNN第二步，查找树中最近的k个点
#
# @param Annoy index, built with AnnoyBuildIndex
# @param query A set of data to be queried against the index
# @param k Number of neighbors
# @param search.k During the query it will inspect up to search_k nodes which
# gives you a run-time tradeoff between better accuracy and speed.
# @param include.distance Include the corresponding distances in the result
#
# @return A list with 'nn.idx' (for each element in 'query', the index of the
# nearest k elements in the index) and 'nn.dists' (the distances of the nearest
# k elements)
#
#' @importFrom future plan
#' @importFrom future.apply future_lapply
#
AnnoySearch <- function(index, query, k, search.k = -1, include.distance = TRUE) {
  #> a$getNTrees()
  # [1] 50 #有50棵树，就是对各个维度的划分方法。

  # (A1)行数: cell 数
  n <- nrow(x = query)

  # (A2)构建空矩阵，n行 k列：保存每个细胞(共n个)的k个最近邻点？
  idx <- matrix(nrow = n,  ncol = k) #放索引？
  dist <- matrix(nrow = n, ncol = k) #放距离

  # (A3) 第一个参数是不是 Rcpp_AnnoyAngular 类的对象? 默认F
  convert <- methods::is(index, "Rcpp_AnnoyAngular")

  # (A4) 如果 plan() 不是 multicore 对象
  if (!inherits(x = plan(), what = "multicore")) {
  	# 则记录现在是单核，退出时恢复这个单核状态。
    oplan <- plan(strategy = "sequential")
    on.exit(plan(oplan), add = TRUE) #退出时执行， add 这个参数干啥的？
  }

  # (A5) 遍历每个细胞  1:n
  # 想要真的多线程，还需要在运行前，手动设置:
  #library(future)
  #library(future.apply)
  #plan(multisession, workers=8) #多线程声明
  #message("outside>>how many cores can use now: ", nbrOfWorkers()) #8

  res <- future_lapply(X = 1:n, FUN = function(x) {
  	# 查找第x=1个细胞，找到最近的k个邻居，
  	# 查找最多 search.k 个节点，默认 -1
  	# 包含距离，默认T
    res <- index$getNNsByVectorList(query[x, ], k, search.k, include.distance)
    
    # Convert from Angular to Cosine distance | 默认是欧拉距离，不需要转换
    if (convert) {
      res$dist <- 0.5 * (res$dist * res$dist)
    }
    
    # 返回一个list，有2个元素
    # C++中下标是0，回到R中+1
    list(res$item + 1, res$distance)
  })

  # (A6) 遍历每个细胞 1:n
  # 对于每一行，填充一个数据：该行最近邻的 k 个行的下标和距离。
  for (i in 1:n) {
  	# 每个list的第一个元素，最近点的编号
    idx[i, ] <- res[[i]][[1]]
    # 每个list的第二个元素，最近点的距离，默认T
    if (include.distance) {
      dist[i, ] <- res[[i]][[2]]
    }
  }
  # (A7) 返回一个list，第一个元素是最近邻的k个点的坐标(n行k列)，第二个元素是到最近邻的k个点的距离
  return(list(nn.idx = idx, nn.dists = dist))
}







()Indices.Neighbor()
./seurat-object-4.0.4/R/generics.R:452:Indices <- function(object, ...) {

#' Get Neighbor nearest neighbor index matrices
#'
#' @param object An object
#' @param ... Arguments passed to other methods;
#'
#' @return A matrix with the nearest neighbor indices
#'
#' @rdname Indices
#' @export Indices
#'
#' @concept data-access
#'
Indices <- function(object, ...) {
  UseMethod(generic = "Indices", object = object)
}




./seurat-object-4.0.4/R/neighbor.R:110:Indices.Neighbor <- function(object, ...) {

#' @rdname Indices
#' @export
#' @method Indices Neighbor
#'
Indices.Neighbor <- function(object, ...) {
  object <- UpdateSlots(object = object) #更新
  #获取每个元素最近邻的k个元素的下标构成的矩阵
  indices <- slot(object = object, name = "nn.idx") 
  # 一行是一个cell，添加cid
  rownames(x = indices) <- slot(object = object, name = "cell.names")
  return(indices)
}







() Graph 的定义 (S4 类)

#' The Graph Class
#'
#' The Graph class inherits from \code{\link[Matrix:sparseMatrix]{dgCMatrix}}.
#' We do this to enable future expandability of graphs.
#'
#' @slot assay.used Optional name of assay used to generate \code{Graph} object
#'
#' @name Graph-class
#' @rdname Graph-class
#' @exportClass Graph
#'
#' @seealso \code{\link[Matrix]{dgCMatrix-class}}
#'
Graph <- setClass(
  Class = 'Graph',
  contains = "dgCMatrix", #定义父类
  slots = list(
    assay.used = 'OptionalCharacter'
  )
)




() as.Graph.matrix()

把稀疏矩阵强转为其自定义子类
按后缀分发方法，这是一个S3泛型函数。

#./seurat-object-4.0.4/R/generics.R:55:as.Graph <- function(x, ...) {
#./seurat-object-4.0.4/R/graph.R:52:as.Graph.Matrix <- function(x, ...) {
#./seurat-object-4.0.4/R/graph.R:75:as.Graph.matrix <- as.Graph.Matrix
#./seurat-object-4.0.4/R/graph.R:84:as.Graph.Neighbor <- function(x, weighted = TRUE, ...) {

as.Graph.Matrix <- function(x, ...) {
  CheckDots(...)
  x <- as.sparse(x = x)
  if (is.null(x = rownames(x = x))) {
    stop("Please provide rownames to the matrix before converting to a Graph.")
  }
  if (is.null(x = colnames(x = x))) {
    stop("Please provide colnames to the matrix before converting to a Graph.")
  }
  return(as(object = x, Class = "Graph"))
}

as.Graph.matrix <- as.Graph.Matrix



测试：
mat <- matrix(data = 1:16, nrow = 4)
rownames(x = mat) <- paste0("feature_", 1:4)
colnames(x = mat) <- paste0("cell_", 1:4)
g=as.Graph(x = mat)


输出:
> mat
          cell_1 cell_2 cell_3 cell_4
feature_1      1      5      9     13
feature_2      2      6     10     14
feature_3      3      7     11     15
feature_4      4      8     12     16

> g
4 x 4 sparse Matrix of class "Graph"
          cell_1 cell_2 cell_3 cell_4
feature_1      1      5      9     13
feature_2      2      6     10     14
feature_3      3      7     11     15
feature_4      4      8     12     16

> class(g)
[1] "Graph"
attr(,"package")
[1] "SeuratObject"
> is(g)
 [1] "Graph"          "dgCMatrix"      "CsparseMatrix"  "dsparseMatrix"  "generalMatrix"  "dCsparseMatrix"
 [7] "AnyMatrix"      "dMatrix"        "sparseMatrix"   "compMatrix"     "Matrix"         "xMatrix"       
[13] "mMatrix"        "replValueSp"   

> g@assay.used
NULL







()ComputeSNN() C++函数
#./seurat-4.1.0/R/clustering.R:968:ComputeSNNwidth <- function(
#./seurat-4.1.0/R/RcppExports.R:96:ComputeSNN <- function(nn_ranked, prune) {

ComputeSNN <- function(nn_ranked, prune) {
    .Call('_Seurat_ComputeSNN', PACKAGE = 'Seurat', nn_ranked, prune)
}

#./seurat-4.1.0/src/RcppExports.cpp:431:    {"_Seurat_ComputeSNN", (DL_FUNC) &_Seurat_ComputeSNN, 2},
{"_Seurat_ComputeSNN", (DL_FUNC) &_Seurat_ComputeSNN, 2},


#./seurat-4.1.0/src/RcppExports.cpp:314:Eigen::SparseMatrix<double> ComputeSNN(Eigen::MatrixXd nn_ranked, double prune);
#./seurat-4.1.0/src/RcppExports.cpp:315:RcppExport SEXP _Seurat_ComputeSNN(SEXP nn_rankedSEXP, SEXP pruneSEXP) {
#./seurat-4.1.0/src/RcppExports.cpp:320:    rcpp_result_gen = Rcpp::wrap(ComputeSNN(nn_ranked, prune));

// ComputeSNN
Eigen::SparseMatrix<double> ComputeSNN(Eigen::MatrixXd nn_ranked, double prune);
RcppExport SEXP _Seurat_ComputeSNN(SEXP nn_rankedSEXP, SEXP pruneSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::traits::input_parameter< Eigen::MatrixXd >::type nn_ranked(nn_rankedSEXP);
    Rcpp::traits::input_parameter< double >::type prune(pruneSEXP);
    rcpp_result_gen = Rcpp::wrap(ComputeSNN(nn_ranked, prune));
    return rcpp_result_gen;
END_RCPP
}


#./seurat-4.1.0/src/snn.h:16:Eigen::SparseMatrix<double> ComputeSNN(Eigen::MatrixXd nn_ranked);
using namespace Rcpp;
Eigen::SparseMatrix<double> ComputeSNN(Eigen::MatrixXd nn_ranked);




#############
# SNN 算法的核心函数，算法解析参考本文的4.2
#############
#./seurat-4.1.0/src/snn.cpp:16:Eigen::SparseMatrix<double> ComputeSNN(Eigen::MatrixXd nn_ranked, double prune) {
// [[Rcpp::export(rng = false)]]
Eigen::SparseMatrix<double> ComputeSNN(Eigen::MatrixXd nn_ranked, double prune) {
  std::vector<T> tripletList; //声明一个向量
  int k = nn_ranked.cols(); //总列数

  tripletList.reserve(nn_ranked.rows() * nn_ranked.cols()); //定义行列，初始化全0的
  for(int j=0; j<nn_ranked.cols(); ++j){ //遍历列，共k个
    for(int i=0; i<nn_ranked.rows(); ++i) { //遍历行，细胞数
      tripletList.push_back(T(i, nn_ranked(i, j) - 1, 1)); //添加列向量: 第i行，第_列(R to C, -1); 该位置值为1
    }
  }

  // 声明稀疏矩阵，行数*行数(细胞数)
  Eigen::SparseMatrix<double> SNN(nn_ranked.rows(), nn_ranked.rows());
  SNN.setFromTriplets(tripletList.begin(), tripletList.end()); //填充0矩阵，一行是一个细胞，一行内其最近邻的k个值1
  // https://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html

  SNN = SNN * (SNN.transpose()); //一个0-1矩阵 乘以 其转置

  // 行数 sm1.outerSize(); // Number of columns (resp. rows) for a column major (resp. row major )
  for (int i=0; i < SNN.outerSize(); ++i){
  	// 内部类遍历该对象的行: https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
    for (Eigen::SparseMatrix<double>::InnerIterator it(SNN, i); it; ++it){
      it.valueRef() = it.value()/(k + (k - it.value()));

      //剪枝，小于这个值 默认 1/15，则设置为0
      if(it.value() < prune){
        it.valueRef() = 0;
      }

    }
  }

  // 稀疏矩阵，减掉0的点。
  SNN.prune(0.0); // actually remove pruned values
  return SNN;
}









() todo 
#./seurat-4.1.0/R/clustering.R:627:FindNeighbors.Assay <- function(
#./seurat-4.1.0/R/clustering.R:671:FindNeighbors.dist <- function(

默认执行不到的函数，先跳过。

